function H = H_jacobian( x,y )
H(1,1)=(x-x_s1)*((x-x_s1)^2+(y-y_s1)^2)^(-1/2);
H(1,2)=(y-y_s1)*((x-x_s1)^2+(y-y_s1)^2)^(-1/2);

H(2,1)=(x-x_s2)*((x-x_s2)^2+(y-y_s2)^2)^(-1/2);
H(2,2)=(y-y_s2)*((x-x_s2)^2+(y-y_s2)^2)^(-1/2);

H(3,1)=(x-x_s3)*((x-x_s3)^2+(y-y_s3)^2)^(-1/2);
H(3,2)=(y-y_s3)*((x-x_s3)^2+(y-y_s3)^2)^(-1/2);

H(4,1)=(x-x_s4)*((x-x_s4)^2+(y-y_s4)^2)^(-1/2);
H(4,2)=(y-y_s4)*((x-x_s4)^2+(y-y_s4)^2)^(-1/2);
end

